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LDDMM SURFACE REGISTRATION WITH ATROPHY CONSTRAINTS 


Abstract. Diffeomorphic registration using optimal control on the diffeomorphism group 
and on shape spaces has become widely used since the development of the Large Defor¬ 
mation Diffeomorphic Metric Mapping (LDDMM) algorithm. More recently, a series 
of algorithms involving sub-riemannian constraints have been introduced, in which the 
velocity fields that control the shapes in the LDDMM framework are constrained in ac¬ 
cordance with a specific deformation model. Here, we extend this setting by considering, 
for the first time, inequality constraints, in order to estimate surface deformations that 
only allow for atrophy, introducing for this purpose an algorithm that uses the augmented 
lagrangian method. We also provide a version of our approach that uses a weaker con¬ 
straint in which only the total volume is forced to decrease. These developments are 
illustrated by numerical experiments on brain data. 


1. Introduction 

Over the last couple of decades, multiple studies have provided evidence of anatomical 
differences between control groups and cognitively impaired groups at the population level, 
for a collection of diseases, including schizophrenia, depression, Huntington’s or dementia 
mmmmmmmmmmm- in the particular case of neuro-degenerative diseases, 
a repeated objective has been to design anatomical biomarkers, measurable from imaging 
data, that would allow for individualized detection and prediction of the disease. This 
goal has become even more relevant with the recent emergence of longitudinal studies, 
involving patients at early stages or “converters” which showed that, when the effect is 
measured at the population level, anatomical changes caused by diseases like Alzheimer’s 
or Huntington’s were happening several years before cognitive impairment can be detected 
on individual subjects. 

Surface registration [21J using the LDDMM algorithm is a powerful tool for the analysis 
of shape variation in ROIs represented by triangulated surfaces. While one its main ad¬ 
vantages is its flexibility and its ability to render smooth, diffeomorphic, free-form shape 
changes, there are situations in which prior information is available, and should be used to 
enhance the results of the shape analysis. In this paper, we focus on situations in which 
no tissue growth is expected to occur (like with brain longitudinal data). In such contexts, 
it is natural to ensure that shape analysis should only detect atrophy, even when noise 
and inaccuracy in the ROI segmentation process may lead in the other direction. (Here, 
we mean “atrophy” in the general sense of local volume loss.) In this paper, we introduce 
an atrophy-constrained registration algorithm, that include some of the ideas introduced in 
[2j, while extending them to inequality constraints associated to the problem we consider. 
This algorithm will be described in section [2j with our numerical approach discussed in 
section [3j Some theoretical results on existence of solutions and consistency of discrete 
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approximations are provided in section [4j An extension of the algorithm to include affine 
alignment and its restriction to scalar volume decrease constraint are provided in section 
[5] and [6j Finally, experimental results are provided in section [7j 

2. Atrophy-constrained LDDMM 

2.1. Continuous Optimal Control Problem. The LDDMM algorithm implements an 
“optimal control” strategy in which a template surface So is “driven” toward a target 
surface Si via a time-dependent process t i-a S(t), with S(0) = So- This is achieved by 
minimizing 

( 1 ) ^I^Wvmldt + DiS^S,) 

subject to the state equation dS/dt = v(t,S(t)) where v is a smooth velocity field on M 3 . 
By HI y, we mean a functional norm in a reproducing kernel Hilbert space (RKHS) V. that 
we will assume to be embedded in Co(R d ,R d ) (the completion, for the standard supremum 
norm of up to p derivatives, of the space of compactly-supported infinitely differentiable 
vector fields) with p> 1. This space can, for example, be defined as the Hilbert completion 
of 

(Av(x)) T v(x)dx 

(originally defined for smooth vector fields), where A is a differential operator. More gen¬ 
erally, letting A : V i —y V* be the Hilbert duality mapping, and K = A -1 , the reproducing 
kernel of V, also denoted K, is a mapping defined on M 3 x M 3 , with values in Ad 3 (H) (the 
set of 3 by 3 real matrices) such that, letting K*(x,y) denote the ith column of K(x,y), 
the vector field K*(-, y) : x rA K?(x. y) belongs to V for all y £ M 3 with, for all v £ V, 
(!€*(•, y ), v) v = Vi(y), the ith coordinate of v(y). To simplify the notation, we will assume 
in this paper that 1C is a scalar kernel, i.e., that it takes the form HC(x,y) = K(x, y) Id]g 3 
where K is-scalar valued. 

The function D in (jTJ) is a measure of discrepancy that penalizes the difference between 
the controlled surface S(-) at the end of its evolution and the target surface S\. Among 
the measures introduced in the literature in combination with the LDDMM algorithm, the 
most convenient computationally are designed as Hilbert space norms between surfaces 
considered as linear forms over spaces of smooth structures. The simplest example, linear 
forms on smooth scalar functions arising from integrating functions over surfaces, yields 
the “measure matching” cost introduced in mm)- “Current matching”, introduced in 
pa m results from integrating smooth differential forms over oriented surfaces. More 
recently varifold-based cost functions j5j were designed, in which functions defined on 
M 3 x Gr(2,M 3 ) (the Grassmannian manifold of 2D spaces in M 3 ) are integrated over the 
surface. Details on these cost functions, their discrete versions on triangulated surfaces, 
and the computation of their gradient are provided in the cited references. 

In optimal control language, v is the “control”, S is the “state”, and v is optimized in 
order to bring the state close to a desired endpoint. With this construction, each point xq 
in Sq is registered to a point x(t) = pit, xq) in S(t) that evolves according to the differential 
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equation dx/dt = v(t,x), with x(0) = xq. The overall evolution is diffeomorphic, i.e., for 
each time t, <p(t , •) can be extended to a smooth invertible transformation with smooth 
inverse on M 3 . 

To define our atrophy constraints, we assume that surfaces are closed and oriented. We 
let No(t,xo) be the outward-pointing unit normal to Sq. An outward-pointing normal to 
5(t) at x = tp(t,x o) is then given by 

(2) N(t, x ) = dip(t, x 0 )~ T N 0 (t, x 0 ) 

where d(p(t,x o) denote the differential of y i-A (p(t,y) at y = xo (a 3 by 3 matrix), with 
the —T exponent indicating the inverse transposed (N(t, •) does not necessarily have norm 
one). 

We will express our atrophy constraint by the fact that the surface evolves inward at all 
points, i.e., by v(t, x) T N(t, x) < 0 for all x E 5(f) and f E [0,1]. Adding this constraint to 
the original surface-matching LDDMM problem leads to the atrophy-constrained problem 

i 

|| v(t) ||ydt + -0(5(1), 5i) -A min 
subject to dtS = v(t,S(t)), v(t, x) T N(t, x) < 0. 



We now reformulate this problem under the assumption that 5o is parametrized with 
an embedding : M -a M 3 , where M is a two-dimensional Riemannian manifold. This 
is no loss of generality, since one can always take M = Sq and go = identity. We take 
parametrizations q : M —> M 3 as state variables, together with functions N : M —> M 3 and 
solve 


Problem 1. Minimize 

9(0,-) = 9o> N (0, •) = A r o, 
d t q(t,-) = v(t,q(t,-)), 
d t N(t, •) = -dv(t, q(t, ■)) r N(t , •), 
v(t,q{t,-)) T N(t,-) < 0 

where, with a slight change of notation, iVo(m) is the outward-pointing unit normal 
to 5o at qo{m). N(t,m) is then an outward-pointing (not necessarily unit) normal to 
5(f) = q(t,M) = ip(t,So) at q(t,m). The third equation in the constraints is the time 
derivative of ([2]). 

2.2. Discrete Approximations. We now assume that surfaces are triangulated, identi¬ 
fying 5 with a pair ( q,F ), where q E (M 3 ) n is a set of n vertices (where n depends on 
5) and F C {l,...,n} 3 lists the indices of vertices that form the triangular faces. We 
assume that the surface is oriented, so that an edge which belongs to two faces is oriented 
in different directions in each face. If q = (q±,..., q n ), and / = ( i,j,k ) E F. the area 
weighted normal to / is 

1 

/) = gfa; ~ ft) x («fc “ ft)> 


(3) i 


\v(t)\\ydt + D(q(l, M), Si), subject to 


(4) 
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which is invariant by circular permutation of i, j and k. From this we define the area- 
weighted normal at vertex q k by 

(5) N k (q,F)= J2 

feF-.kef 

with N(q, F) = (Ni(q, F ),..., N n {q, F)). 

To define a discrete version of Problem [l] we introduce a small relaxation parameter 
e > 0 and state variables q = (qi,..., q n ) and N = (Ah,..., N n ), initialized with a initial 
surface So = (q 0 ,Fo), with a target surface Si = (q 1 ,Fi). The discrete problem will 
minimize 


1 

2 



v(t)\\ydt + Si), subject to 


'q( 0, •) = q 0 ,N(0) = N{q 0 ,F 0 ), 
dtqk(t) = v(t,q k (t)), 
d t N k (t ) = -dv(t,q k (t)) T N k (t), 

^(t^kit)) ■ N k (t ) < e|AT fe (t)|, k = l,...,n 


with S(t) = (q(t), Fq). The introduction of the parameter e is justihed by Theorem [2j and 
is also natural given discrete approximation errors. 


Given that the end-point cost and the constraints depends on v only through the con¬ 
figurations q(t), one shows, using standard RKHS reductions, that the optimal v takes the 
form 

n 

v(t,-) = ^2 K (-,q k (t))a k (t) 

k= 1 

where K is the reproducing kernel of V. Using this, the previous problem reduces to 

Problem 2. Minimize 

i r i n 

(6) - ^ K(q k (t),qi(t))a k (t) ■ ai(t)dt +D(S(l),Si) 

" / ° k,l =1 


(7) 


subject to 


' q{ 0 , •) = q 0 ,N(0) = N{q 0 ,F 0 ), 

n 

d t q k (t) = ^2 K ® (*))«/(*)> 

i=i 

n 

9tN k (t ) = ~^2diK(q k (t),qi(t))N k (t) ■ a t (t), 

i=i 

-Ez=i K iqk(t),qi{t))ai(t) ■ N k (t) < e\N k (t)\,k = 1,... ,n 


with S(t) = (q(t),Fo). 


However, in the discrete case, it is possible to avoid the introduction of A?" as a state 
variable and solve instead: 
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Problem 3. Minimize 


( 8 ) 

(9) 


subject to 


i r l n 
■'° fc,i=i 


' 9(0,-) = g 0; 

n 

dt Qk(t) = ^2 K (q k {t),qi(t))ai(t), 

l=i 

n 

^ ~2K(q k (t),qi(t))ai(t ) ■ N k (q(t),F 0 ) < e|IV fc (<?(*),F 0 )|> 
„ z=l 


with S(t) = ( q(t),Fo ). 


k = 1,..., n 


Note that the apparent simplification is balanced by the fact that the constraint becomes 
a more complex function of the state and controls, with N k (q, Fq) given by 0 and ([5]). 


3. Numerical Method 

Problem [3] is solved using augmented Lagrangian methods, introducing, as described in 
US], slack variables to complete inequality constraints. More precisely, let 

C k i{q) = K(q k ,qi)N k (q,F 0 ) T 

and C(q) = ( C k i(q ), k, l = 1,..., n) the associated n x 3 n matrix. Let K(q) be the 3n x 3n 
matrix formed from blocks ( K(q k , qi) Id K 3, k,l = 1,..., n). For a vector u, let u + denote 
the vector formed with the positive parts of each of the coordinates of u. Let |JV| denote 
the n-dimensional vector will kth coordinates equal to \N k \. 

The augmented Lagrangian is defined by 


(10) F(a, A) 



a(t) T K (q(t))a(t)dt + D(S(1),M) 



(c{q(t))a{t) 


e\N(t)\ 



2 

dt 



|A(f)| 2 df. 


Here, the parameter /r is a small positive number, A E M n is a Lagrange multiplier and 
q is considered as a function of cc via the state equation dtq = K(q)a. The augmented 
Lagrangian optimization procedure iterates the following steps (starting with initial values 
(c*o, Ao)): 

(1) Minimize a 1 — > F(a,X n ) to obtain a new value a n+ i (and q n _ n via the state 
evolution equation). 

(2) Update A, with A n+ i = -// (C(g n+1 )a re+ i - e\N\ - ^iA„) + . 

(3) If needed (e.g., if | (C(q n+ i)a n+ i — e|Af|) + | is larger than a threshold 6 n ), replace 
H by a smaller number, pn, with p < 1. 
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The most expensive step is, of course, the first one, which requires to solve an optimal 
control problem similar in complexity to the unconstrained problem. The computation of 
\7 a F(a, A) (the gradient of F with respect to a) uses a back-propagation algorithm, also 
called the adjoint method. Since similar computations are described in multiple places 

1231123 El El, we briefly summarize it here. 

(1) Given a, compute the associated state q via dtq = K(q)a and the matrix C(q). 

(2) Introduce a co-state p(t ) = (pi,... ,p n ) G (M 3 )", t G [0,1], defined by 
p(l) = V<j.D((g(l), Tq), Si) and 


dtp = -V 9 | p T K(q)a - J a T K(q)a - | 


C(q)a — er|iV|- 


(3) Define V a F = K(q)(a - p) + p [( C(q)a. - e\N\ - — ) ) C(q). 


4. Existence of Solutions and Convergence 

Theorem 1. Assume that V is continuously embedded in Cq(M 3 ,M 3 ) for p > 2 and 
that the data attachment term D is such that <p i-A D(<p(So), Si), defined over all C p - 
diffeomorphisms ip such that p—id£ C f , is continuous for the uniform C p convergence over 
compact sets. Then Problems 01 andjs] always have an optimal solution v G L 2 ( 0,1, V). 

Proof. The theorem is proved along the same lines as similar statements addressing the 
existence of solutions for LDDMM problems [201 3 GEL 23 G] • Considering Problem [lj 
and starting with a minimizing sequence v n , the boundedness of v n in L 2 ( 0,1; V) implies 
that (extracting a subsequence if needed) v n converges weakly in this Hilbert space to a 
limit v (with a norm in L 2 ( 0,1; V) no larger than the liminf of the norms of v n ). This, 
in turn, implies that the associated flows of diffeomorphisms p n (associated to v n ) and 
their first p spatial derivatives, converge uniformly (in time and space) over compact sets 
to (p (associated to v) and its first p derivatives. One gets from this that D(p n ( 1, So), Si) 
converges to D(ip(So), Si). Letting ( q n ,N n ) denote the state in Problem 1 associated to 
the control v n , we have q n (t) = p n (t) o qo —>• q. This implies that the cost function is 
minimized at v. Note that N n (t) = dp n (t)~ T Nq and therefore converges to N(t). 

We now show that v(t,q(t,x)) T N(t,x) < 0 for all x G M and almost all t G [0,1]. 
Fixing x, one easily deduces from the facts that v n converges weakly to v, dv n is uniformly 
bounded, q n and N n converge uniformly to q and N, that v n (t, q n (t , x)) T N n (t, x) converge 
weakly to v(t, q(t, x)) T N(t , x) in L 2 (0,1; M), which implies that v(t, q(t , x)) T N(t, x) < 0 for 
almost all t (the set of non-positive a.e functions is a closed convex set in L 2 and therefore 
also weakly convex). By considering a countable dense set of x's, and using the fact that 
x i —> v(t, q(t, x)) T N(t, x) is continuous, we find that v(t, q(t, x)) T N(t, x) < 0 for all x G M 
and almost all t G [0,1]. The same argument can be used to prove existence of solutions 
for Problems [2] and [3 □ 
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Next, we adress the problem of convergence as the triangulation is refined. More pre¬ 
cisely, as the triangulation of a fixed initial surface So gets finer and finer, do solutions of 
Problems [2] and [3] converge to solutions of Problem [lj We will see that we need to slightly 
relax the constraints. 


Theorem 2. Assume that V is continuously embedded in C^M 3 ,!! 3 ) for p > 2, with 
IMIv < cmax. rgR 3(|u(x)|, \dv(x)\) where |u(.t)| is the Euclidean norm of v(x) and \dv(x)\ 
is the operator norm of dv(x) and c > 0 is fixed. Also assume that if a sequence of 
triangulations (S n ) ne ^ converges to a surface S in the sense of currents, then D(S n ,S i) —> 
D(S,S i). 

Let (So)neN = (qo,T n ) n > 3 be an increasingly fine triangulation with n points of a com¬ 
pact surface So of class C 2 with normal vector field N. We assume q k q = q k 2 0 for any 
k < min(ni,n 2 ) (so that q k , k E N can be defined independently of n), that the averaged 
normals Nf 0 at q k0 converge to N(q k0 ) after normalization as n goes to infinity, and that 
{qk, k £ N} is dense in So- Then, there exists a decreasing sequence or real numbers e n > 0 
with e n —>■ 0, such that if v n £ L 2 ( 0,1; V) solves Problem [i] with e = e n , then the sequence 
(u n ) ne N is bounded and any weak limit point v* of v n is a solution to Problem [7} 


Proof. First of all, since ||u n (f)|| 2 < 2D(Sq, Si) and D(Sq , Si) —> D(So, Si), we see that 
(v n ) n £f$ is bounded in L 2 { 0,1; V ) and therefore contained in a weak compact subset. We 
can assume that (u n ) ne pj is weakly convergent without loss of generality. Let v* be the 
weak limit of ( v n ). We start by proving that v* satisfies the constraints of Problem [lj 
Let ( q n ,N n ) denote the states for the discrete problem. Also let q(.{t) satisfy dtq(.{t) = 
v*(t, qk{t)) with <7j*(0) = go,fc> dtS*(t ) = v*(t , S*(t )) with S*(0) = S°, and N*(t ) be the unit 
normal vector field to S*(t). 

Then, using the same method as in the previous theorem, we get that dtqft = v(-, qj((-)) 
converges weakly to dtqjl in L 2 and that Nf(-)/\Nf(-)\ converges strongly to N*(t,q(.(t )) 
in C°. As a direct consequence, for every / : [0,1] —> M smooth nonnegative function, we 
easily obtain d t q* k {t) ■ N*(t , q* k (t))f(t)dt < 0, so that d t q* k {t) ■ N*(t, q* k {t )) = v*(t, q%(t)) ■ 
N*(t,q k (t )) < 0 for almost every t and every k. Therefore, v* satisfies the constraints of 
Problem [U 

We still need to prove that v* is optimal. Let v be a solution of Problem [lj Then t i—>- 
||u(£)|| is constant and no greater than ^2D(So, Si). Let <p be the flow of v. Using Gron- 

dif~ T (t,q n )N n , 

wall’s estimates, one obtains a sequence e n —> 0 such that v(t, ip(t , qQ k )) ■ \d^~ T {t )N^\ — 

e n , with £ n depending only on V, D(Sq,Si), and the value of N(q k 0 ) — Nff 0 /\N k 0 
Hence, v satisfies the constraints of the discrete problems, so that 


1 

2 



v n (t)\\ 2 ydt + D(S n (l),S 1 ) < 



v(t)\\ydt + D(ip(1,Sq),Si). 


Letting n go to infinity, we get the same inequality with v* and S* in place of v n and S n , 
so that v* is optimal. 
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The value of e n can be explicitely bounded from above with respect to the C 2 -norm of 
a parametrization of Sq. Note that this proof can be adapted to the case of Problem [3j 
although the exact values of e n may differ. □ 


5. Affine Alignment 

Because the RKHS V is embedded in Cq(M 3 ,M 3 ), vector fields v G V vanish at infinity. 
This implies, in particular, that affine transformations do not belong to this Hilbert space, 
and that the diffeomorphic registration does not incorporate any rigid alignment. If such 
an alignment is needed, one can include it in the optimal control framework by completing 
the control with the corresponding vector fields. 

Let the registration be computed over GkI 3 , where G is a subgroup of GL-j, (R), ix 
referring to the semi-direct product extending G with translations to obtain affine trans¬ 
formations. Let g be the Lie algebra of G, with basis E\...., Eh. Instead of v G V, we use 
a control given by (v, /3i,..., (3h, r) with (3\,... ,/3h 6 l and r G M 3 , and the state equation 

h 

( 11 ) d t q k {t) = v(t,q(t)) + J £2f3i(t)Eiq(t) + r(t) 

i=i 


with associated cost \ J' ( j f||v(t)||y + J2k=l c kPk(t) 2 + c ol r (^)| 2 J dt for some non-negative 
numbers Co, c \,..., c^. The extension of the numerical procedure to this setting is rather 
straightforward, and not detailed here. Of course, the affine components must be added to 
v in the atrophy constraint v ■ N < 0. 

Consider the special case G = SO 3 is the rotation group (so that h = 3) and assume 
that Euclidean transformations act as isometries on V, which means that, for all v G V, 
R G SO 3 , b G M 3 , the vector field v : x i-a R t v(Rx + b) also belongs to V and ||h||y = ||v||y. 
Euclidean-invariant RKHS’s of vector fields are extensively described in |Sj], to which we 
refer for more details. In the case of scalar kernels K(x, y) = K(x,y)Id^ 3 , Euclidean 
invariance implies that K is a radial kernel, i.e., that I\(x,y) = ^(\x — y | 2 ) for some 
function 7 (additional conditions on 7 are needed to ensure that K is a positive kernel; see 
0 ). Assume finally that the end-point cost D is invariant by Euclidean transformation: 
D(S, S') = D(RS + b, RS' + b). In this case, the optimal control problem (using cq = ■ ■ ■ = 
C 3 = 0 ) that minimizes \ f ( ] ||u(f)||ydt + D(S( 1), Si) in v, (3, r, subject to 

q( 0 ) = qo, 

d t q{t) = v(t, q(t)) + X)?=i Pl{t)Eiq{t) + r(t) 
d t N(t, ■) = - (dv(t, q(t, ■)) + E?=i Pi(t)Ei) N(t, ■) 

(v{t, q(t)) + Ef=l Pl{t)Eiq(t) + r(t )) • N(q(t)) < 0 


(12) 
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is equivalent to minimizing 2 J Q \\v(t)\\ydt + D(S(1), RiS\ + b\) in v, R\, b±, subject to 

q{ 0 ) = Qo , 


(13) 


d t q(t) = v(t,q(t)) 
d t N(t, •) = -dv(t,q(t,-)) T N(t, 
K v(t,q(t)) ■ N(t,q(t )) < 0 


via the change of variable q(t) = R(t)q(t ) + 6 (f), v(t,x) = f?(t) 1 v{t,Rx + 6 ), IV(f) = 
R(t)N(t), with d t R = Ef =1 d tb = E?=i + r{t), Ri = f?(l ) _1 and 

b\ = —i?(l)^ 1 6 (l). In other terms, Euclidean registration via optimal control using (12) is 


equivalent to the original atrophy-constrained LDDMM optimizing its target over an orbit 
under the action of rigid transformations. 

Note that co = ■ ■ ■ = Q, = 0 should not be used with general affine transformations 
when non-compact components of the affine group are added to rotations. Intuitively, 
this would allow small deformations to be scaled up to larger ones at zero cost, and one 
can conjecture that the associated optimal control problem has no solutions, and admits 
minimizing sequences of controls with vanishing cost at the limit. The equivalence with 
a problem in which the target is allowed to vary over its orbit via affine transformations 
is not true either for groups larger than the Euclidean one, essentially because invariant 
kernels do not exist in such cases. 


Some attention should be paid to the time discretization, in particular in the rigid case. 
In our implementation, we use a simple Euler scheme to discretize the equation dtq = v(t, q), 
i.e., we take q(t + 5t) = q(t) +5t.v(t , q(t)). When using rigid registration, however, we take, 
with A(t) = X)f=i Pl{t)Ei a skew-symmetric matrix 

q(t + 5t ) = e St q(t) + 5tv{t , q(t )) + 5tr(t) 

to discretize dtq = v(t, q)+Aq+r, with the explicit expression e L> = Id + S1 ” Cu U + 1 ~ c ° sCu U 2 , 

° u c u 

c u = y/—tr(U 2 ) for a 3 by 3 skew-symmetric matrix U. This ensures that the rigid 
registration remains a displacement, even for large values of the coefficients /?/, which are 
made possible by the absence of cost on this transformation. 

6. Global Volume Constraint 

Ensuring that the total volume of the surface decreases (instead of enforcing inward 
motion at every point) can be done very similarly to the full atrophy constraint, using the 
single constraint J2k =l v ( t ^ % W)' -Njfe(i) < 0 or, after reduction Y%,1 =l K (qk(t),qi{t))ai(t) ■ 
-Nfc(f) < 0, where Nf. is given by ([5]). It is important here to use the area-weighted normal, 
to discretize the continuous constraint J 50 ) v (t, ■) ■ A < e, which provides the 
derivative of the total volume with respect to time (where volg( t ) i s the volume form 
on 5(f)). This constraint can be rewritten as l^C(q(t))cst(t) < e, where l n is the n- 
dimensional vector with all coordinates equal to 1. The reformulation of the augmented 
Lagrangian method for this scalar constraint is straightforward and left to the reader. 
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7. Experimental Results 

Fig. [l] provides three examples of segmented hippocampus surfaces taken from the BIO¬ 
CARD longitudinal study H2EU- Each row corresponds to a different subject, starts with 
a baseline image, and then compares the outputs of unconstrained, volume-constrained and 
atrophy-constrained surface registration when mapping to a follow-up surface. The color 
map is proportional to the total normal displacement during the estimated deformation. 
The segmented follow-up of the first two subjects is slightly larger in volume than the 
baseline and this is corrected by the two constrained method. The deformation pattern 
in the volume-constrained approach is however very similar to the unconstrained case (the 
deformed surface is only slightly, and almost uniformly, smaller in the constrained case, 
no ensure that no volume is gained). The deformation pattern in the atrophy-constrained 
case is completely different, since it can only include atrophy (blue). The last subject had 
significant volume lost, so that the volume-constrained and the unconstrained registrations 
return identical results, both having a few regions with outward growth, which disappear 
in the atrophy-constrained case. 



Figure 1. Follow-up image is shown in wireframe in all views. From left 
to right: baseline; LDDMM surface registration (no constraint); registration 
with non-increasing volume constraint; registration with atrophy constraint 
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